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Morphological properties of strained epitaxial films are examined through a mesoscopic approach 
, developed to incorporate both the film crystalline structure and standard continuum theory. Film 

04 ' surface profiles and properties, such as surface energy, liquid-solid miscibility gap and interface 

}_( ' thickness, are determined as a function of misfit strains and film elastic modulus. We analyze the 

I stress-driven instability of film surface morphology that leads to the formation of strained islands. 

We find a universal scaling relationship between the island size and misfit strain which shows a 
crossover from the well-known continuum elasticity result at the weak strain to a behavior governed 
I by a "perfect" lattice relajcation condition. The strain at which the crossover occurs is shown to be 

a function of liquid-solid interfacial thickness, and an asymmetry between tensile and compressive 
strains is observed. The film instability is found to be accompanied by mode coupling of the complex 
amplitudes of the surface morphological profile, a factor associated with the crystalline nature of 
O . the strained film but absent in conventional continuum theory. 

I. INTRODUCTION 

The most recent area of focus in thin film epitaxy has been on exploiting the growth and control of strained solid 
[ films to develop specific nanostructure features that can be used in optoelectronic device applications. These structures 
include junctions, quantum wells, and multilayers/superlattices for which planar interfaces are highly desired. On 
the other hand, epitaxially grown films are usually strained due to the lattice mismatch with the substrate, leading 
to a variety of stress-induced effects and structures either on the film surface or across the interfaces, such as islands 
Q (quantum dots) or nanowires.^^'* A wide range of device applications results from such heterostructures, including 
Q , LEDs, diode lasers, detectors, FETs, etc.,^''^ with the major technical concerns being the requirement of long-range 
ordering, size regularity, placement and defect control. 

Much progress has been made in understanding film growth above the surface roughening temperature, particu- 
larly the formation and evolution of coherent nanostructures. The evolution sequence often involves many physical 
processes, including an initial morphological instability of the Asaro-Tiller-Grinfeld (ATG) type^"— that results in 
surface ripples and undulations , ^ ^ d^ the formation of islands and the evolution from pre-pyramid to faceted shape 
' (e.g., {105}-faceted pyramids for SiGe^^), subsequent islands coarsening, 1^"- further shape transitions from pyramids 
to domesi^ or to unfaceted prepyramids^^ and the nucleation of misfit dislocations for very large islands 
fT^ ■ To understand these complex processes of nanostructure self-assembly, most of current theoretical efforts are based 
I on either continuum diffusion and elasticity theories or atomistic simulation methods that focus on a certain single scale 
of description. In standard continuum theory, the film morphology is described by a coarse-grained, continuum surface 
profile*'^ or phase fields, ^^^■^^ with evolution governed by the relaxation of continuum elastic and surface free energies. 
Quantitative results have been obtained to reveal fundamental mechanisms of film nanostructure formation observed 
in a variety of experimental systems. Recent work has focused on morphological instabilities of strained films^i"— 
^ , or superlatticesr^"— the coupling to alloy film composition inhomogeneity^"^ island evolutionj^Si^ ordering and 
' coarsening^ ^^^ as well as island growth on nanomembranes/nanoribbonsj '^^i'^^ Such continuum approaches give 
a long- wavelength description of the system, which has a large computational advantage over microscopic approaches 
but naturally neglects many microscopic crystalline details that can have a significant impact on film structural 
evolution and defect dynamics. This can be remedied via atomistic simulations such as kinetic Monte Carlo (MC) 
methods. Recent progress includes identifying detailed properties of strained islands such as morphology, density and 
size distributioii^i^^ and the evolution of complex surface structures including dots, pits and grooves as a function of 
growth conditions in both two^ and three^ dimensions. However to simulate strained film growth, novel approaches 
(e.g., Green's function metho d'^^'^^ or local approximation technique'*^) are required to incorporate strain energy via 
long-range elastic interactions, which usually limit atomistic studies to small length and time scales. 

Recently an approach coined Phase Field Crystal (PEG) modeling has been developed to incorporate atomic-level 
crystalline structures into standard continuum theory for pure and binary systems.—^— This model can be related to 
other continuum field theories such as classical density functional theory^^— and the atomic density function theoryi^ 
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The PFC model describes the diffusive, large-time-scale dynamics of the atomic number density field p, which is spa- 
tially periodic on atomic length scales. By including atomic scale variations, the physics associated with elasticity, 
plasticity, multiple crystal orientations and anisotropic properties (of, e.g., surface energy and elastic constants) is nat- 
urally incorporated. This approach has been applied to a wide variety of phenomena including glass formation^ climb 
and glide dynamics of dislocations)^ epitaxial growth^^^ii^ii^i^^i^ pre-melting at grain boundaries j^i^ commensu- 
rate/incommensurate transitionS )^^i^° sliding friction phenomena'^ and the yield strength of polycrystals.— i^>^i^ 
For strained film epitaxy, the basic sequence of film evolution observed in experiments, i.e., morphological instability 
— > nanostructure/island formation — > dislocation nucleation and climb, has been successfully reproduced in PFC 
simulations Unfortunately computational simulations of the original PFC model are limited by the need 
to resolve atomic length scales. This limitation can be overcome by deriving the corresponding amplitude equation 
formalism as developed by Goldenfeld et al^^ to effectively describe the system via the "slow" -scale amplitude and 
phase of the atomic density p, while at the same time retaining the key characters (e.g., elasticity, plasticity and mul- 
tiple crystal orientations) of the modeling. Very recently such a mesoscopic approach has been extended by Yeon et 
al^ to incorporate a slowly-varying average density field which is essential to account for the liquid-solid coexistence 
and a miscibility gap, and also by Elder et alM- to describe the binary alloy systems for both two-dimensional (2D) 
hexagonal and three-dimensional (3D) bcc and fee structures. Application of this extended expansion to strained film 
growth and island formation has yielded promising results, particularly the determination of a universal size scaling 
of surface nanostructurcs (strained islands).^** However, in these PFC studies some key factors for understanding the 
basic mechanisms of strained film evolution are still missing and yet to be addressed, including film surface properties 
(such as strain-dependent surface tension and width) and the effect of the sign of film/substrate misfit strain, as will 
be clarified in this work. 

In this paper we provide a complete formulation for such multiple-scale analysis of single-component, strained film 
epitaxy. Compared to our previous worl^ which is also based on the amplitude equation formalism established for 
two-dimensional high temperature growth, here we provide a new and more systematic study of various strained film 
properties including surface energy, film surface (or liquid- film interface) thickness, and liquid-film miscibility gap 
that are identified for different misfit strains (both tensile and compressive). Furthermore, morphological instabilities 
of the strained films and the corresponding behavior of island formation are systematically investigated, showing the 
important effects of misfit strains (both magnitude and sign) and film surface properties that are absent in previous 
work. A main feature of our multi-scale (mesoscopic/microscopic) approach is that it can maintain the efficiency 
advantage of the continuum theory through coarse-grained amplitudes, without losing significant effects due to the 
discrete nature of the crystalline film structure. 



II. AMPLITUDE EQUATION FORMALISM FOR STRAINED FILM EPITAXY 

In the PFC model;^'^'^ the free energy functional F can be derived from the classical density functional theory 
of freezing*^ and be expressed in terms of a dimensionless atomic number density n = {p — p)/p, i.e., 

F/pkBT ^ J dr{^[B' + B- (2i?2v2 + i?4y4)] n ~ ^n' + ^n^} , (1) 

where p is the average density, T is the temperature, R represents the lattice spacing, B^ is related to the isothermal 
compressibility of the liquid phase, B^ is proportional to the bulk modulus of the crystalline state, and r and v 
are phenomenological parameters (chosen as r = 1/2, v — 1/3 in the following calculations for simplicity). The 
liquid-solid transition is controlled by a parameter e = {B^ — B^)/B^ which is related to temperature difference from 
the melting point. The solid phase exists at e > 0, with hexagonal/triangular crystalline symmetry in 2D and bcc in 
3D. Based on the assumption of conserved system dynamics, i.e., dn/dt — TW^SF/Sn with P the mobility, the PFC 
dynamic equation is given by 

dn/dt = PV^ [B'^n + ^^(i?^^' + 2R^\/^)n - rn^ + vn^] . (2) 

Defining a length scale Iq — R, a. time scale tq — R^ /TB^ , and n — > ^Jv/B^ n, we obtain the rescaled equation 

dn/dt = [-e n + (V^ + ql)^n - gn^ + n^] , (3) 

where g = r /\J vB^, qo — 1 and the symbol qo is retained for the clarity of presentation. 

For the epitaxial system of interest, we consider a system configuration composed of a semi- infinite strained crys- 
talline film and a coexisting homogeneous liquid state, which are separated by a time-evolving interface (i.e., film 
surface) . To access the "slow" time and length scales of the film surface profile we introduce a standard multiple scale 
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expansion of the PFC equation ([3]) and derive the associated amphtude equations, with detailed procedures given in 
Refs. l66l - [68l For a 2D system with the film surface normal to the y direction, the atomic density field n is expanded 
in both liquid and solid regions as the superposition of a spatially /temporally- varying average local density no (for 
the zero wavenumber mode) and three hexagonal base modes, i.e., 

3 

n = no(^, Y,T) + Y^ Aj{X, Y, T)e<-'- + c.c, (4) 

where both no and complex amplitudes Aj are slowly varying variables (with Aj = in the liquid region), and 
represent the three hexagonal basic wave vectors 

/ \/3. ~ /^V3~ ... 

Qi = 90 I — ~ 2^ I ' " ^° 1 "2"^ " 2^ I ■ 

This expansion (U) implies the separation of "slow" scales X = e^/'^x, Y = e^^'^j/, T — et for no and Aj (and hence 
the film surface profile) from the underlying crystalline structure, at the limit of small e or high temperature growth. 
The corresponding amplitude equations are given by (in the form of Model G^) 

dAj/dt = -ql6T/6A*, (6) 
dno/dt = V^6T/6no, (7) 

where the effective potential J- (a Lyapunov functional) is written as 

33 „ 3 



J- = fdr\i-e + 3nl - 25no) ^ | + ^ | (V^ + 2zq° . V) + ^ ^ \A 
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Note that the operator (V^ +2iq"- V) preserves the rotational covariance of these amplitude equations. This effective 
free energy describes a first order phase transition from a liquid (Ai = 0) to a solid state {Aj ^ 0) and incorporates 
elasticity though the operator (V^ + 2zq^ • V), as discussed in Ref. [691 In addition the terms containing no incorporate 
a miscibility gap for the density at liquid-solid coexistence. 

For a hexagonal lattice, the equilibrium wave numbers along x and y directions are q^^ = V3qo/2 and qy„ — go for 
the undistorted, zero- misfit bulk lattice. For strained films during epitaxy (with distorted hexagons/triangles), the 
misfit Sm is determined by 

e,„ = ^ = ^-l, (9) 
a qx„ 

where = 27: /qx^ is the stress-free bulk film lattice constant and a = 27: /q^ is the lattice constant of the strained 
film. The complex amplitudes Aj should then be expressed by 

where amplitudes A'^ are complex, 5^ = Qxo^m — \/iqo£m/2, and the value of 6y (7^ S^) is determined by the lattice 
relaxation along the film growth direction y (corresponding to the Poisson relaxation in continuum elasticity theory) . 
Since both Aj and A', are slowly varying quantities, Sx, Sy and the misfit strain (e^) should also be sufficiently small. 
Substituting Eq. (fTOl) into Eqs. (jG])-®, the amplitude equations for strained films are then 

dtA[ = -ql { [-e + 3n^ - 2.gno + {d^ + 3^ - i{V3qo + 2Sx)dx - i{qo + Sy)dy 

-VSqodx -Si- qody/2 - dl/^y A[ + (6no - 2g)A'*Ai,* 

+3A[{\A[\^ + 2\A',\^ + 2\A',\^)}, (11) 



dtA'2 = -ql { [-e + 3n2 - 25no + {dl + + 2z((jo + Sy)dy - 2qoSy - S. 
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These amplitude equations describe a strained system and will be used to study morphological instabilities of a 
liquid-crystal surface under strain. In the next section, steady state or base solutions will be obtained for a planar 
liquid-crystal interface under strain. In Sec. IIVI the stability of these planar solutions to small perturbations at the 
surface will be examined. 



III. BASE STATE SOLUTION: FILM SURFACE PROPERTIES 



We first construct a base state involving a planar film surface (i.e., a coexisting liquid-crystal interface). The 
corresponding amplitudes A^, and density nj] are then only a function of the normal direction y, and hence the 
amplitude equations ([TT |) - p3)) can be simplified as 



dAydt = ~ql5F^/6Af, dnl/dt ^ d^SJ^ySn] 



(15) 
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dy - i{qo + Sy)dy - VSqoS., -si- qoSy/2 - 52/4] Al 

[d^y + 2l{qa + 5y)dy - 2q„6y - 6^] A^]^ 
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(6n° - 29){AlAlAl + Af) + 6 {\A',mf ' 
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The equilibrium profile for the base state (with solid/liquid coexistence) is given in Fig. [TJ corresponding to non- 
growing, stationary films of different misfit strains £m and elastic constants (as determined by B'-^). The amplitudes 
and Tip can be used to reconstruct the full density field n via Eq. ([4]), as shown in Fig. [2l This figure highlights the 
increase in interfacial width as the magnitude of elastic moduli (i.e., B^) increases. Since the stationary solution of 
Eqs. ()15p and (|16p cannot be obtained analytically, the results shown were obtained by numerical solutions based 
on a pseudospectral method. To apply the periodic boundary condition, we set the initial configuration as a pair 
of symmetric liquid-solid interfaces located sX y — Ly/4, and ZLy/A respectively, with Ly the one-dimensional (ID) 
system size which is chosen up to Ly = 8192 in our calculations so that these two interfaces are sufficiently far apart 
from each other and thus evolve independently. In the numerical algorithm adopted, the second order Crank- Nicholson 
time stepping scheme is used for the linear terms, while a second order Adams-Bashford explicit method is applied for 
the nonlinearities. A grid spacing Ay — Ao/8 (i.e., 8 grid points per basic wavelength Aq = 2TT/qo) is chosen in most 
of calculations, although similar results have been obtained with much larger Ay. Relatively large time steps At can 
be adopted without losing numerical stability: We use At = 0.5 (or even 1) for > 10, and At = 0.2 for i?^ — 1 
with sharp interface. We also use the same algorithm and parameters in the stability /perturbation calculations given 
in Sec. |IV1 

For finite misfits the amplitudes \Ai \ = jAgl ^ jvljl and their difference increases with Em as shown in Fig. [H This 
corresponds to a triangular structure distorted along the y direction (the surface normal) and the degree of distortion 
increases with misfit strain. Also as shown in Fig. [TJ for larger value of B^ which corresponds to smaller bulk modulus 
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FIG. 1: The equilibrium (solid/liquid coexistence) profile of the base state, for e = 0.02, — 1 and 10, and misfit Em = 
(solid fines) and 5% (dashed lines). 
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FIG. 2: Sample equilibrium profiles of the complete density field n as reconstructed from Uq and Aj for e — 0.02, a misfit of 
3%, and = 1 and 10 in a) and b) respectively. 



(as we calculate based on one- mode approximation; see Sec. lIVp . the interface or film surface is more diffuse (i.e., 
with larger interface width), but with a narrower coexistence region (i.e., smaller but nonzero miscibility gap). This 
can also be seen in Fig. |4l which shows the liquidus and solidus rescaled density well as the miscibility 

gap AriQ — nQ°' — n^'^ as a function of misfit Em- The size of miscibility gap decreases with the increasing magnitude 
of misfit, and shows slight asymmetry with respect to the misfit sign as a result of different non-linear elastic effects 
on liquid-solid coexistence property for tensile and compressive strains. 

We also calculate the surface tension 7 as a function of misfit strain since it is one of the important factors for 
determining film stability and island formation. Surface energy is known to play a stabilization role on film evolution 
and for simplicity is often approximated as misfit independent in many strained film studies 9- 2^ However in the 
presence of a strain field, the surface energy is known to vary as a result of intrinsic surface stress tr" and is usually 
expanded up to 2nd order in terms of strain tensor Uij (with i, j the film surface coordinate indices) in linear elasticity 
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FIG. 4: (a) The equilibrium densities n^'^ and no°' (in the coexisting liquid and solid regions respectively) as a function of 
misfit strain Sm, with parameters the same as those of Fig. [3] (b) The size of miscibility gap Ano = no°' — rig'^ as a function 
of Sm, for = 10 (the main panel) and 1 (the inset); Note the large vertical scale for = 1 in the inset. 



theory,2iZ^ i.e., 

7 = 7o + <^ijUij + ^SijkiUijUki, (17) 

where Sijki are the surface excess elastic moduli. Both tT^^- and Sijki can be either positive or negative.— For the 
ID surface considered here, strain u^x = Em and hence Eq. (|T7| gives 7 = 70 + cf^x^m + Sxxxxs'inl'^i which is 
consistent with our amplitude-equation calculations shown in Fig. [S] Data fitting of our numerical results yields 
70 = 6.82 X 10-3, cr°^ = -4.77 x 10-^ Sxxxx/2 = -9.76 x IQ-^ for = 1, and 70 = 2.20 x IQ-^, = -3.72 x IQ-^, 
Sxxxx/"^ = —2.06 x IQ-^ for i?^ — 10 (all in dimensionless unit), showing smaller surface energy for larger value of 
(with larger surface width). These results indicate that for the parameters chosen, both the intrinsic surface stress 
CT°^ and excess elastic moduli Sxxxx are negative, leading to the decrease of surface energy with increasing magnitude 
of misfit strain. In addition the tensile surface stress is rather weak which can explain the weak asymmetry of 7 
between tensile and compressive strained films. 
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FIG. 5: Results of surface tension 7 for different misfit strains Em, with parameters the same as those of Fig. [3l Also shown 
are the quadratic fitting results for = 10 (in the main panel) and 1 (in the inset). Note that the vertical scale in the inset 
for = 1 is much larger. 



IV. MORPHOLOGICAL INSTABILITY AND ISLAND SCALING 



For strained films with nonzero misfit, a morphological instability of film surface is known to develop as a result 
of strain energy relaxation, leading to surface undulations and then the formation of surface nanostructures such as 
strained islands. Such an instability can be revealed via a linear analysis of amplitude equations given above. We can 
expand the amplitudes in Fourier series as 



no{x,y,t) = n'^iy) + J2no{q^,y,t)e'-'^-'^ , 

<3i 



(18) 
(19) 



where A'^(y) and riQ^y) are the planar base solutions discussed in the previous section and the perturbed quantities 
Aj and ho obey the following linearized equations, 

dtAi{qx,y, t) = -ql { + inf - 2gnl + (d^ - i{qo + 5y)dy - ql + (TSgo + 2(5x)gx 
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The stability of the base planar film surface is examined by introducing initial small random perturbations into 
Aj and no, and solving numerically the initial value problem defined by Eqs. (|20p - ((23)) . given a specific value of qx- 
The numerical algorithm introduced in Sec. IIIII is employed, with the use of a pseudospectral method and periodic 
boundary conditions. 

For nonzero misfit, within a certain range of wave number qx the initial perturbations of Aj and no grow with 
time around the liquid-solid interface, while they always decay to zero far from the interface region, showing the 
stability of both the solid and liquid bulks. This interface instability results in the formation of islands or mounds 
at the liquid-solid interface, as shown in Fig. [51 This figure was obtained by reconstructing full density field n 
from the amplitudes with wave number qx of maximum instability (based on Eq. (HI). A typical example of the 
dynamics of the amplitudes that gives rise to this instability is given in Fig. [7^. We then calculate the perturbation 
growth rate cr{qx), noting that |no| e*^*. This process is repeated for a range of perturbation wave number 
qx, and also for various misfits Em- Some results of the dispersion relation are shown in Fig. [7)3, for e = 0.02 and 
= 10. Previous work of continuum elasticity or phase-field theory has predicted various forms of dispersion relation, 
including a ~ aag^ — a4g^ (for surface-diffusion dominated processr'"— ) a ~ —a2q^ + ct-^q^ — a^q^ (if considering 
wetting effectsr^i2^) a ~ a\q — a-i'f' (in the case of evaporation-condensationjSiiSiiSS) or ct ~ a-iq^ — 035^ (for bulk- 
diffusion dominated case^^) with q the wave number and ai {i — 1, ...,4) the model-dependent coefficients that are 
usually a function of surface tension and elastic moduli. However, none of these forms fits our dispersion data, which 
instead can be well fitted only by a 4th order polynomial of qx for all range of wave numbers, similar to a combination 
of all the above forms. This is not unexpected, given that all factors of surface diffusion, bulk diffusion, wetting 
effects, and evaporation/condensation are naturally incorporated in the PFC model and cannot be easily decoupled. 
This can be seen through the fact that the PFC modeling of epitaxial growth involves the coexistence of liquid-solid 
interface that buckles and evolves, and thus naturally involves the diffusion processes along the interface and between 
liquid region and solid film, and also the variation of material properties such as surface/interface energy and elastic 
relaxation across the interface (i.e., the wetting effects). We expect that an important parameter controlling these 
different processes would be e, the temperature distance from the melting point. The e (or temperature) dependence 
of properties of system relaxation has been known for pattern formation systems, and is also seen in our PFC studies. 
Here we focus on high temperature regime where the amplitude equation representation is most relevant and effective, 
and hence choose e = 0.02 which is different from other studies with larger e and hence lower growth temperature 
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FIG. 6: Reconstruction of full density field n for an interface profile showing island formation, with a 3% misfit at e = 0.02. 
a) corresponds to density n at t = 125, 000, for = 1 and the maximum instability wave number Qx — 0.0184, while b) 
corresponds to n at f = 2000, for B"^ = 10 and = 0.026. 

(e.g., e = 0.1 in Ref. [IH)- For such small e (high temperature) surface diffusion process is more prominent and coupled 
with the bulk diffusion process, a phenomenon that might be weakened or absent in low temperature growth (e.g., 
in Ref. [s^ only bulk diffusion behavior has been identified in the dispersion relation obtained from the original PFC 
equation) . 

The development of surface perturbations and instability can be characterized by an evolution time scale t, which 
can be approximated via the inverse of maximum perturbation growth rate (Xmax and is found to scale as or 
in continuum elasticity theory with the assumed mass transport mechanism dominated by surface diffusion or 
evaporation-condensation respectively.®'^*^ However, our calculations yield results more complicated than this single 
power law behavior, as shown in Fig. [71;, which can also be expected from the coupling of various mass transport 
processes in this modeling as discussed above. Our results show that the time scale r decreases with misfit strain 
Em, since the £„j provides the driving force for the morphological instability, r is also found to significantly decreases 
when increases. For example at a given misfit, r is typically one or two orders of magnitude larger for 5^ = 1 
compared with = 10. This difference is most likely due to the significant decrease in surface energy and increase 
in interfacial thickness as B^ is increased, as shown in Fig. [5] and Fig. [1] respectively. 

The maximum of the growth rate determines the characteristic wave number Qi for the instability, and hence 
the characteristic wave number of the island/mound formation on the film surface. We plot in Fig. [5^ the relation 
of this instability/island wave number Qj vs. misfit strain £„i, for different values of B^ and for both compressive 
{sm > 0) and tensile (e™ < 0) films. For each value of B^ we can identify two regions, corresponding to a quadratic 
behavior of Qi ~ at small misfits (see also the inset of Fig. [SJd) and a linear dependence of Qi on Em for large 
enough strains. Such quadratic scaling in the small misfit limit is consistent with the well-known results of continuum 
theory including all different assumptions of dominant mechanisms such as surface diffusion, evaporation-condensation 
and wetting effect Si^^—i^i^'2^ However, this scaling result differs from the experimental findings in SiGe/Si(001) 
growthfiiii^ which indicate the linear behavior Qj ~ Sm for the stress-driven surface instability and coherent epitaxial 
islands. Although this observation of a linear relationship is qualitatively similar to what we obtain above for large 
enough misfits, it should be cautioned that the experimental systems involve more complicated factors related to the 
SiGe alloying nature that is not considered here, particularly the atomic mobility difference between the two film 
components which was verified by recent first principle calculationsZi and was believed to play a key role on island 
size scaling i^iS 

For the single-component films studied here the crossover from the quadratic scaling at the continuum weak-strain 
limit to linear behavior at high strains is most likely due to the discrete nature of the crystalline lattice that is implicitly 
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FIG. 7: (a) Amplitude perturbations, which grow with time around the solid/Hquid interface for e — 0.02, = 10, wave 
number qx — 0.026 and 3% misfit, (b) Perturbation growth rate a as a function of wave number q^, for different values of 
misfit Em- Other parameters are the same as (a), (c) Characteristic time scale r (~ l/(Tmax) for the mounding instability as a 
function of misfit Em, for = 1 and 10. Two power laws, r ~ and ~ are also shown for comparison. 



included in the amplitude formulation. It is known (and verified in direct simulations of PFC Eq. ([g]) 43,45, 54 -^ ^j^g^^ 
at late times the instability to form islands or mounds leads to the nucleation of dislocations around the edges of 
islands or in the valleys between the mounds. These dislocations nucleate to relieve strain in the film and appear 
at earlier times for larger misfit strains. Here we define a length scale, A^, for "perfect" relaxation such that if the 
dislocations nucleate at this distance apart, strain in the film will be completely relieved (aside from the strain induced 
by the dislocations themselves). We can then make the assumption that if the continuum prediction for most unstable 
wavelength is smaller than A/j, continuum theory will break down. To evaluate Xr consider a 1+1 dimensional film; 
assume being the lateral length of film surface and by definition we have = Na = Muq, where N is the number 
of atoms in strained lattice, M is the atom number for unstrained state after dislocations nucleate, and a and oq are 
the corresponding lattice constants already defined in Eq. Thus from Eq. ^ for the definition of misfit, we obtain 
£m = {N — M)/M, leading to the average distance between dislocations Xr = L/\N — M| = L/{M\em\) = ao/\£m\, 
with the associated wave number Qr — gxo|em| (plotted as a dashed line in Fig. |5^). Assuming that on average 
at least one dislocation will appear at each island edge/valley, this wave number Qr will then be the upper limit 
imposed by the discrete nature of the lattice, as it would be unphysical for islands with size smaller than Xr to appear 
which would instead cause the "overrelaxation" of the film lattice. Our results of island wave number Qj for different 
values of (= 10, 20, 100) all converge to this limit at large misfit strains (except for = 1 which will be discussed 
below). 

This "perfect" relaxation condition is expected to be met at large enough misfits, but not at small strains where 
dislocations appear at far late stage after islands form, leading to the crossover phenomenon between two scaling 
regimes given in Fig. |S1 This crossover occurs when Qi{oi small misfit limit) = Qr. As stated above, at small Em we 



11 




FIG. 8: (a) Characteristic wave number Qi of film surface instability as a function of misfit strain magnitude |em |, for different 
values of = 1 and 10 and both compressive (sm > 0) and tensile (e™ < 0) films. The limit imposed by "perfect" relaxation 
condition is indicated by a dashed line, (b) Scaling of island wave number based on a crossover wave number Q* = 87^0 /4i5 
and misfit = 37^0/4-/?, for different values of which is proportional to film elastic modulus. The inset highlights the 
crossover to the continuum result of Qi ~ at small misfit limit. 



can recover the result of continuum theory which predicts Qj cx {E/j)e'^ (with E the Young's modulus) i^^— In our 
calculations based on the PFC model and amplitude equations, we evaluate E from a one-mode approximation ,^'^i^^ 
E = where Anin = 4(g— 3rto + -\/5^ + 24no5 — 36no + 15£m)/15. Using the results of 7 given in Sec. Hill we 

can fit the small misfit data well into a form Qi = AEe'^/^j (for all values of i?^; see the inset of Fig. |5}d). Therefore, 
the misfit (ej^) and island wave number {Q*j) at the crossover can be determined via Q*j = AEe^^ /i^ — Qn = q^o^mj 
resulting in = i^q^o/AE and Q\ — i^q^Q/AE. Defining rescaled quantities Q — Qi/Q*i and im — ^m/Sm^ we can 
then scale all the data from different conditions (e.g., films of different elastic constants, for > 1) onto a single 
universal scaling curve accommodating all range of misfit strains, for both compressive and tensile films (see Fig. 
IHb). The crossover misfit strain can be very small (< 2%, depending on e.g., film elastic properties), showing the 
breakdown of continuum approach even at relatively large scales. 

Note that although this linear behavior due to "perfect" lattice relaxation and the scaling crossover have been 
observed in our previous work,— it was limited to compressive strained films and not-too-large misfits. However, the 
more generalized study given here shows a small deviation from the limit of "perfect" relaxation for small value of 
, as indicated in Fig. [8^ with island wave numbers of B^ = 1 lying above such upper limit (the dashed line) when 
the magnitude of mismatch \em\ exceeds 5% (for tensile films) or 6% (compressive). Similar deviation can be seen in 
the corresponding scaling plot of Fig. ^jp. Nevertheless, at large misfits the linear scaling behavior is still maintained, 
which is qualitatively different from the quadratic scaling at the small strain limit. Based on the discussions given 
above for "perfect" relaxation condition, it is expected that Qi > Qr occurs only when some of the island edges 
would be dislocation-free even at late evolution times. The condition for this scenario is not clear; but our results 
suggest that this may occur when the liquid-solid interface (or film surface) is sharp enough. As given in Fig. [HI 
the interface width W decreases with the value of i?^, and is particularly small at = 1 (with W ^ 13.5Aj/ for 
both tensile and compressive films, less than 2 lattice spacing) as compared to others. It could then be expected that 
details of film morphological evolution, including instability and island formation, would be different for such sharp 
interface, as somewhat indicated in Fig. |S1 Further studies are needed to clarify this special scenario of strained film 
evolution. 

Fig. [S] also yields the effect of finite interface width W on the island size (or wave number) scaling. We find 
I/Q7 i.e., a linear relation between crossover instability wavelength (= 2ti/Q*j) and the interface thickness. 

This is consistent with most recent results of direct PFC simulations^^- which indicate that the discrepancy or crossover 
between the classical elasticity result of quadratic scaling of Qi and the linear behavior identified in the PFC modeling 
could be attributed to the finite thickness of the interface, a fact that is neglected in the classical continuum theory. 
As seen in Fig. [HI when — >■ (i.e., the assumption adopted in continuum elasticity theory), we have Q*j ^ 00 
and hence recover the continuum theory prediction of Qi ^ for the whole range of misfit strain, as expected. 
Corresponding to real experimental systems. Fig. [9] predicts that at constant growth temperature (same e value), the 
liquid-solid interface thickness varies with film elastic modulus (or the value of -B^), and for different film materials 
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FIG. 9: The linear relation between the liquid-solid interface width W and the inverse of crossover wave number 1/Q/, for 
e = 0.02 and both compressive (filled symbols) and tensile (open symbols) films. Values of are also indicated on the plot. 



the crossover island size separating two island scaling regimes increases linearly with the interface thickness. 

Another important feature of our results is the asymmetry between tensile and compressive films which, however, 
becomes distinct only at small enough and large enough misfits (see Fig. [8] for the data of = 1). Given the 
important role played by the surface energy 7 on film stability and evolution, we expect this asymmetric phenomenon 
of island wave number to be closely related to the property of 7 shown in Fig. [S] The intrinsic surface stress a^^ 
determined for B^ = 1 is an order of magnitude larger than that for B^ — 10, leading to much larger value of surface 
energy difference between tensile and compressive strains; also such difference increases with the magnitude of misfit 
strain. The corresponding behavior of surface instability and island formation would then follow the similar trend, as 
observed in Fig. [51 



V. FREE ENERGY ANALYSIS AND MODE COUPLING 



To further elucidate the properties of the strained surface, it is interesting to analyze the effective free energy J- 
(given in Eq. ([5])). Consider the net change of J- relative to that of a planar interface, i.e.. 



(24) 



where is the free energy of the planar interface given in Eq. ([TBI) . AT < indicates film surface instability against 
the initial perturbation, while AJ^ > refers to the energy penalty of any perturbations and hence corresponds to 
stability of planar film surface. 

Based on the Fourier expansion and , AT can be expanded up to second order in the perturbed quantities 
Aj and ho, i.e., 



(25) 



Detailed expression of the first order term AJ^^^' is given in the Appendix (see Eq. ([Al[) ). We find numerically 
AJ-^^^ ^ 0, and hence the net energy change AJ- is determined by the second order quantity 



AJ-(2) ^AT+ + AT_, 



(26) 



where 



A^_ =L,Jdyj:\ (6nO - 2g) ^ ^°i*(-fc) + AfA,{q,) hl{q,) 
+ {&nl - 2g) A°Ai{q^)A2{-qa:) + A°Ai{<lx)A3{-qa:) + A°A2{'lx)A3{-q^) + c.c. } 



(27) 
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FIG. 10: Time evolution of effective free energy density change AJ^ (per unit volume V — L^Ly) of the perturbed state, with 
misfit Em ~ 3%, wave number Qi = 0.026, e = 0.02, and -B^ = 10. Also included are the positive contribution AJ^+ and the 
negative contribution AT-. 



with Aj{qx) = Aj{qx,y,t) and no(9£c) = naiQx,y,t), and the contribution AJ^^ is shown in Eq. (IA2p of the Appendix. 

Given the numerical solution for the perturbed amplitudes (see Eqs. (pOjl -ip ^ ') as described in Sec. IIVI AJ^ 
(~ AJ^*^^^) can be approximated via the most unstable characteristic wave number by substituting the numerical 
solutions for amplitudes at Qx — ±Qi- We find that all terms in Eq. (IA2p are positive, i.e., AT+ > 0; both two 
terms in Eq. (P7| yield negative contribution (noting that usually Qjiq ~ 2g < for liquid-solid coexistence), so that 
AJ^_ < 0, and the magnitude of the last term is much larger than the 1st one. As shown in Fig. [TUl at large enough 
time A dominates over the stabilizing terms in AJ-^ , leading to negative net free energy change AJ- and thus the 
film instability. Note that the last term in Eq. (I27p , which dominates AJ-- , arises from the 2nd-order expansion of 
{A1A2A3 + A^A^A^) in the effective free energy formula ([5]). It represents the coupling of different modes of complex 
amplitudes, and our numerical results show that it contributes to the integral of AJ^_ only in the interface or film 
surface region (as the perturbed amplitudes decay fast in the bulks). We can then argue that it is the mode coupling 
of complex amplitudes at the liquid-solid interface that is mainly responsible for the morphological instability of the 
strained film. Note that the amplitudes of structural profile Aj are complex, and thus their evolution involves an 
important process of phase perturbation (or phase winding) . Physically this phase behavior corresponds to the elastic 
relaxation of the lattice structure, and thus the mode coupling property identified above indicates that the coupling 
of elastic relaxation for different lattice modes (or wave vectors) around the film surface would be one of the major 
factors underlying the film instability and mounding behavior. Such phase behavior is related to details of crystalline 
structure, as captured by the PFC model and the amplitude equation formalism, but not by the continuum theory. 
Furthermore, the competition between AJ^_|_ (> 0) and AJ^_ (< 0) shown in Fig. [10] is consistent with previous 
analysis of continuum elasticity theory showing the competition between film stabilization effects (such as surface 
energy) and destabilizing factors (mainly elastic effects) 4"— i^"— Note also that the above mechanism identified 
should be already incorporated in the original PFC equation ([3]) and the associated PFC free energy ([T|), while the 
analysis given here based on the amplitude formulation has the advantage of being able to single out individual 
contributions from different lattice modes. 



VI. CONCLUSIONS 



We have investigated the detailed properties of a strained film surface, its morphological instability, and the asso- 
ciated island wave number scaling through a systematic analysis of the amplitude equation formalism based on the 
phase field crystal model. We identify the amplitude and average density profiles of liquid-film coexisting interface, 
the interface width, miscibility gap, and surface energy (including intrinsic surface stress and excess elastic modulus), 
for various misfit strains (both magnitude and sign) and film elastic constants (or values of B^). The morphological 
or mounding instability of the strained film is systematically examined, showing results absent in all previous contin- 
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uum elasticity and phase-field approaches and atomistic modeling. In particular, we obtain a crossover phenomenon 
of instability or island wave number scaling, from the well-known continuum, ATG result of Qj to a linear 

behavior Qj ~ at large enough strains which is identified by an upper limit imposed by the condition of "perfect" 
lattice relaxation. Most data (of different parameter ranges) can be scaled onto a universal scaling relation for the 
whole range of misfit strain, with some small deviations for very narrow liquid-solid interfaces in the large strain 
limit. The asymmetry of film properties between tensile and compressive strains is also observed. Note that although 
either linear or quadratic scaling has been reported in experiments (such as SiGe/Si(001)) and model simulations 
(e.g., kinetic MC) or continuum theory (e.g., ATG instability), the universal scaling relation with crossover of the two 
regions has not been found before. We expect our prediction here to be examined by experiments of single-component 
film epitaxy or atomistic simulations with large enough length and time scales. 

Our study highlights an important feature of the amplitude formulation for strained film epitaxy, in that it can 
simultaneously reproduce continuum results (e.g., the ATG instability) and reveal significant corrections due to the 
microscopic nature of the crystalline structure. Our approach adopts a mesoscopic-level description of the system, via 
the amplitudes or envelopes of the slowly varying surface profile for which the well-developed continuum, mesoscopic 
theory can be applied. On the other hand, the crystalline nature of the strained film is preserved particularly via phase 
perturbations of the complex amplitudes that are prominent around the film surface. The latter has been emphasized 
through revealing the breakdown of traditional continuum approaches even at relatively small misfit stress and the 
associated crossover effect of island size scaling, and also through examining the origin of film instability that is 
accompanied by mode coupling of complex amplitudes in the liquid-solid interface region. Our results thus emphasize 
the importance of multiple scale modeling of complex material systems such as the strained film epitaxy process 
studied above. Note that although in this paper we focus on 2D hexagonal/triangular crystal structure, we expect 
the approach and analysis technique developed here to be directly extended for other crystalline symmetries and 
other surface directions, such as the epitaxial growth and island formation in 3D bcc or fee films for which we have 
developed the corresponding amplitude expansion formulation very recently;^ 
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In this appendix the detailed expansion 
AJ^f-^' shown in Eq. (p5)) . we have 



Appendix A: Free energy expansion 

forms of free energy difference AJ- are presented. For the first order term 



A-F(i) ^L^Jdylj2{-' + 3"o' " 2gn^ + 3|A°p) (^f i,(0) 



+ c.c 



3 



+(6n°-25)^|A°pno(0) 
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+6 + (Ari2(0) + C.C.) + 6 + \Al\^) (^r^(O) + cc.) 



(Al) 



with Aj{0) = Aj{qx = 0,y,t) and no(0) = rio(fe = 0,y,t). For the second order terms, the contribution AJ^+ is given 

by 



+ I [^y - «(9o + 5y)dy -ql + {VSqo + 26x)qx 
-VSqoS, -Si- qo6y/2 - S^/a] 

[dl + 2i{qo + Sy)dy -ql- 2qoSy - si] i2(gx) 
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